Cite this Article
Li Ming-Ru, Zhang Nan, Zhang Feng-Shou. A simulation study of water property changes using geometrical alteration in SPC/E. Chinese Physics B, 2018, 27(8): 083103  
Permissions
A simulation study of water property changes using geometrical alteration in SPC/E
Li Ming-Ru1, 2, Zhang Nan1, 2, Zhang Feng-Shou1, 2, 3, †
Key Laboratory of Beam Technology and Materials Modification of the Ministry of Education of China and College of Nuclear Science and Technology, Beijing Normal University, Beijing 100875, China
Beijing Radiation Center, Beijing 100875, China
Center of Theoretical Nuclear Physics, National Laboratory of Heavy Ion Accelerator of Lanzhou, Lanzhou 730000, China

 

† Corresponding author. E-mail: fszhang@bnu.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11635003, 11025524, and 11161130520), the National Basic Research Program of China (Grant No. 2010CB832903), and the European Commission’s 7th Framework Programme (Fp7-PEOPLE-2010-IRSES) (Grant Agreement Project No. 269131).

Abstract

We present a systematic investigation of the impact of changing the geometry structure of the SPC/E water model by performing a series of molecular dynamic simulations at 1 bar (1 bar = 105 Pa) and 298.15 K. The geometric modification includes altering the H–O–H angle range from 90° to 115° and modifying the O–H length range from 0.90 Å to 1.10 Å in the SPC/E model. The former is achieved by keeping the dipole moment constant by modifying the O–H length, while in the latter only the O–H length is changed. With the larger bond length and angle, we find that the liquid shows a strong quadrupole interaction and high tetrahedral structure order parameter, resulting in the enhancement of the network structure of the liquid. When the bond length or angle is reduced, the hydrogen bond lifetime and self-diffusion constant decrease due to the weakening of the intermolecular interaction. We find that modifying the water molecular bond length leading to the variation of the intermolecular interaction strength is more intensive than changing the bond angle. Through calculating the average reduced density gradient and thermal fluctuation index, it is found that the scope of vdW interaction with neighbouring water molecules is inversely proportional to the change of the bond length and angle. The effect is mainly due to a significant change of the hydrogen bond network. To study the effect of water models as a solvent whose geometry has been modified, the solutions of ions in different solvent environments are examined by introducing NaCl. During the dissolving process, NaCl ions are ideally dissolved in SPC/E water and bond with natural water more easily than with other solvent models.

PACS: 31.15.xv;31.70.-f;61.25.Em
Keyword:molecular dynamics simulation;geometric modification;intermolecular interaction;ion’s solvation
1. Introduction

Water is the most important liquid for our existence. Despite the tiny size of the water molecule, a network structure can be made by strong intermolecular hydrogen bonds (H-bonds) in the liquid ambient condition.[1–4] This capacity to form H-bonds, except as the non-directional interactions seen in some simple liquids, results in many unusual properties. For instance, density exhibits as maximum at 4°C and viscosity shows decreases under high pressure.[5–8] These property changes are strongly related to the structure of the water liquid phase, which in turn depend on the intermolecular potential. Additionally, in aqueous solutions, water as a solvent also has a significant effect on the conformation of proteins and nucleic acids molecule, such as pleated sheets and helices.[9–12] Despite significant research into properties of water, its structural evolution is still somewhat unclear.[13–15] The problem is, how does a water molecule evolve to an optimal geometric structure? How does the intermolecular interaction change and is the solubility of water changed as the water molecule evolves to the optimal structure? Therefore, to explore the structural evolution of the water molecule is very important and intriguing. With the development of computer technology, molecular dynamics (MD) simulations become an important tool: they are not limited by the constraints of experiments, thus offering an ideal route for studying the individual contributions of different aspects.

Two types of modified water models used in such investigations include the alteration of either the geometric angle (Bent model) or the relative strength of electronic and dispersive forces alteration (Hybrid model). Lynden-Bell and colleagues performed a series of molecular dynamic simulations, which resulted in various interpretations regarding the properties of these modified liquids.[16–21] They did so by changing the SPC/E model’s H–O–H angles from 90° to 60° and adjusting the Lennard-Jones (LJ) parameters σ to produce approximately the same dipole moment parameters as the SPC/E water model. This resulted in the water’s network structure changing from a three-dimensional to a one-dimensional chain. The H-bonding strength was affected by altering the weight of the Lennard-Jones term relative to the electrostatic term. This leads to a change of the liquid structure from water-like to a new type of Lennard-Jones liquid without an evident local order structure.

Yet another study[22] focused on investigating the effects of altering the SPC/E model geometric angle as a function of temperature within the angle range from 120° to 60°. They found that the density anomalies appear as the bond angle is at least 100°. Increasing the bond angle leads to a shift of the density maximum anomaly. In addition, their study also showed that the first-neighbour shell for 60° is relatively changed over 200, 250, and 300 K; however, both 90° and 100° gradually acquire water-like structural features upon cooling.

It is inevitable that the change of a water molecule’s structure, and the hydrogen and oxygen atom charge lead to changing its properties. Nevertheless, the dissolving of charged ions strongly depends on the structure or property of the solvent around them. In our pervious research,[23] the relative weight of Lennard-Jones and electronic components was modified, based on the SPC water model, to investigate the relationship between the solution of NaCl and the H-bonding ability of water at room temperature and standard atmospheric pressure. The results indicated that NaCl is most ideally dissolved in normal water with the strongest hydration effects around both cations and anions.

The current study seeks to probe the property changes of water when the geometry structure of the water molecule is altered, and also tries to understand changes in the intermolecular interaction between water molecules during the geometric evolution of a water molecule. So we build a series of water models (not real water) based on the SPC/E water model. Similarly, the solvation behavior of the Na+ and Cl− ions in these water models is also investigated. The O–H lengths variation ranges from 0.90 Å to 1.10 Å (1 Å in SPC/E), while the change in H–O–H angles ranges from 90° to 115° (109.47° in SPC/E). The optimal distance between atoms is found to vary with the change in the SPC/E model geometric structure. The distance change, however, has a significant effect on the strength of H-bonding due to the variation of the inter-molecular interaction. The average reduced density gradient (aRDG) and thermal fluctuation index (TFI) methods are used further to study the change of intermolecular interaction in a dynamic environment and provide a good way to distinguish H-bond, vdW interaction, and steric effect. Finally, we consider the solvation of NaCl in these models. These studies can help us to further gain fundamental understanding of the microscopic origin of the water molecule during geometric evolution.

The remainder of this paper is organized as follows. In Section 2, we outline the simulation procedure used. We summarize and discuss the simulation results in Section 3. Finally, our conclusions are given in Section 4.

2. Model and methods
2.1. Model

The SPC/E water model[24] is used for the simulation in conjunction with classical molecular dynamics simulations. The force field parameters of the SPC/E model, Na+ and Cl− ion are presented in Table 1. The dispersive Lennard–Jones interaction of the SPC/E water model is only provided by its oxygen atom. The geometric modification models that the atomic charge and the Lennard-Jones parameters are kept consistent in the original SPC/E model. All parameters of the geometric modification models are summarized in Table 2. In our models, the model (I) that the molecule’s dipole moment keeps the same 2.35 Debye by changing the O–H distance from 0.82 Å (90°) to 1.08 Å (115°). For the models (II), however, the O–H length is merely altered because the molecule’s dipole moment changes slightly from 2.12 Debye (0.90 Å) to 2.59 Debye (1.10 Å). In what follows, the Bent models and Elastic models are described by simplified expressions, e.g., the bond angle 115° model is denoted by B115, and the bond length 1.10 Å model by E1.10.

Table 1.

Values of Lennard-Jones and electrostatic interaction potential parameters. e represents the magnitude of the electronic charge.

.
Table 2.

Modified water models. ∠HOH, rOH, and rHH are respectively the water molecule’s bond angle, bond length, and distance between two hydrogen atoms. Model (I) represents the change of the H–O–H angle. Model (II) represents the change of the O–H length.

.

The potential energy E between two molecules is written as a sum of pairwise intermolecular interactions and can be expressed as

where Elj describes the Lennard-Jones interaction between two molecules, and Ec is the electrostatic interaction between two charged sites.

2.2. Simulation details

In our simulations, the following molecular dynamic simulation protocol is adopted. Every simulation cell contains 256 water molecules for pure solvent in a cubic box. The v-rescale thermostat[26] and Berendsen barostat[27] are used to regulate the temperature T at 298.15 K and pressure P at 1 atm (1 atm = 1.01325 × 105 Pa), with time constants of 0.2 ps and 0.5 ps, respectively. The leapfrog algorithm with a time step size 1 fs and the SHAKE algorithm[28] are used respectively to integrate the equation of motion and constrain the intramolecular geometry for the rigid models. The cutoff radius of the Lennard-Jones and Coulombic interactions is 9 Å. The Particle Mesh Ewald (PME) method is used to treat electrostatic interactions.[29,30] In the ion solution, the molarity of NaCl is 2.3 M. Every simulation cell contains 234 solvent molecules, 11 Na+, and 11 Cl− ions with cubic periodic boundary conditions. To simulate the variations introduced by ions, we use the same ensemble conditions as pure solvents. Finally, the resulting 70-ns NPT simulation is carried out for each case using GROMACS software.[31]

2.3. Reduced density gradient (RDG) and thermal fluctuation index (TFI)

The reduced density gradient (RDG) and the thermal fluctuation index (TFI) are carried out by the Multiwfn program.[32] The RDG can be written as follows:

where ρ(r) is the electron density and |∇ρ(r)| is the electron density gradient mode. The RDG is a dimensionless quantity describing the difference between the actual electron density and a homogeneous electron distribution. For the average reduced density gradient (aRDG), the average electron density and average electron density gradient mode can be obtained through calculating the dynamics trajectory.

Furthermore, the stability of weak interaction between molecules can be given by the TFI based on the aRDG and is defined as:

where std [ρ(r)] is the standard deviation of electron density in a dynamical trajectory, which can be written as:
where n is the number of frames in consideration and ρi(r) is the density calculated based on the geometry of frame i. All the different models iso-surfaces are plotted with the VMD code.[33]

2.4. Hydrogen bond lifetimes

The H-bond lifetimes can be computed in different ways.[34,35] In the present study, we use a geometrical criterion to determine the existence of H-bonding as shown in Fig. 1. When rOO ≤ 3.5 Å and α ≤ 30°, the water pair is considered as an H-bond.

Fig. 1. Geometrical hydrogen bond criterion.

The value of 3.5 Å corresponds to the first minimum of the radial distribution functions (RDFs) of the SPC/E water model. The lifetime of the H-bonds is calculated from the average over all autocorrelation functions of the existence functions (either 0 or 1) of all H-bonds:

with si(t) = {0, 1} for hydrogen bond i at time t. The integral of C(τ) gives an estimate of the average H-bonding lifetime τHB:

3. Results and discussion
3.1. Structure

The local structure of liquid can be conveniently studied by the radial distribution functions (RDFs). In Fig. 2 and Fig. 3, a set of RDFs for various water models is plotted, while in Fig. 4 a series of the oxygen spatial distribution functions (SDFs) from the different models gives complementary information to the RDFs.

Fig. 2. The RDFs of (a) oxygen–oxygen and (b) oxygen–hydrogen atom pairs for various bond length models.
Fig. 3. The RDFs of (a) oxygen–oxygen and (b) oxygen–hydrogen atom pairs for various bond angle models.
Fig. 4. (color online) Spatial distribution functions of local oxygen density around the water molecule in the first shell in the various water models. It shows an isovalue of 10.5 for B115, compared to an isovalue of 5.5 for other water models.

As shown in Fig. 2 and Fig. 3, the enhanced peak and the deep valleys are exhibited by a larger H–O–H angle and a longer O–H length. This implies that the network structure of liquid becomes more structured than in the SPC/E model. Furthermore, the distance between atoms of two adjacent water molecules can be changed by altering the geometrical structure, resulting in a strong effect on the intermolecular interaction. Specifically, we can see the shift with respect to the first peak in Fig. 2(b) and Fig. 3(b), the optimal O···H distance between oxygen–hydrogen atoms are increased with decreasing O–H length or H–O–H angle, hence leading to the strength of the H-bonding being declined. Figure 2(b) shows the shift of O···H distance transforms from 2.0 Å (E0.90) to 1.6 Å (E1.10), while Figure 3(b) shows the O···H distance shifting from 2.2 Å (B90) to 1.6 Å (B115).]. Similarly, the optimal O···O distances of two adjacent water molecules are also increased with decreasing the O–H length or H–O–H angle, resulting in the weakening of the effective inter-oxygen attraction. Figure 2(a) shows the O···O distance shift from 2.9 Å (E90) to 2.6 Å (E110), while figure 3(a) indicates a shift from 2.9 Å (B90) to 2.7 Å (B115). Significantly, as for Fig. 2(a) and Fig. 3(a), as the bond length and angle decrease, their second peaks gradually fall out. Finally, the value of the peak is twice that of the first shell at E0.90 and B90, which is a character of Lennard-Jones liquid.

Svishchev and Kusalik[36] used firstly the SDFs to study the local structure of water. Figure 4 shows their results for the structures of the first shell around a water molecule for different water models. Especially, the isovalue is the ratio of local molecules density to average density and the central molecule has been included to define the local frame.

As seen in Fig. 4, although all the models exhibit a tetrahedral arrangement, the variation of the ratio of local molecules density to average density is noticeable in bulk water. This variation is caused by the change of the geometry structure of water. For the B115 and E1.04 models, both models indicate an exact tetrahedral structure arrangement with neighboring water molecules. This result further points to the network structure being enhanced when the O–H length and H–O–H angle are increased. Additionally, figure 4 clearly exhibits the variation of the different models as a consequence of decreasing the H–O–H angle and O–H length, where the area behind the oxygen atom the local density is gradually declined. The effect is probably mainly due to a significant reduction of the intermolecular interaction strength leading to a decrease in attraction between water molecules.

3.2. Local tetrahedral order and quadrupole interaction

The parameters of orientational order ⟨q⟩ have often been used in recent years to investigate structural order in liquids and glasses.[37–40] The is defined by

where ψjk is the angle between the bond vectors rij and rik, with j and k representing the four nearest tetrahedral atoms, such that perfect tetrahedra corresponds to q = 1. In addition, these model’s quadrupole interactions QT can be measured using[41,42]
where 2θ is the H–O–H angle and μ is the dipole moment.

As shown in Fig. 5, QT increases with increasing ⟨q⟩ in both types of geometrically modified water models. QT is proportional to the strength of the quadrupolar interactions. For B115 and E1.10 models, the tetrahedral structure order parameter increases to 0.81 and 0.84, the corresponding values of QT are 2.55 Debye · Å and 2.46 Debye · Å, respectively [SPC/E’s ⟨q ⟩ and QT are respectively 0.63 Debye · Å and 2.35 Debye · Å]. This implies that both of the liquids form an ice crystal structure due to the enhancement of the intermolecular interaction.[37,43]

Fig. 5. Variation of the local order parameters ⟨q⟩ and quadrupolar interactions QT for various H–O–H angles, and O–H length models.

The bond angle or length decrease, due to the weakening of the intermolecular attraction, resulting in the water molecule becoming more mobile. This indicates that the local tetrahedral order of water molecules is disrupted. Thus, ⟨q⟩ is reduced monotonically. In addition to the variation in gOO(r) [shown in Fig. 2(a) and Fig. 3(a)] and ⟨q⟩, it can be seen that the liquid changes from a high tetrahedrally ordered state to an LJ-like liquid as the geometric modification decreases monotonically.

3.3. Self-diffusion constants and hydrogen-bond lifetime

The self-diffusion constants of various models were calculated from their mean squared displacements

with the correlation times ranging from 0 to 1 ns.

The correlation of self-diffusion constants and H-bonding lifetime for various H–O–H angles and O–H length models are depicted in Fig. 6, clearly showing that the longer the H-bonding lifetime corresponds to the smaller the self-diffusion constant. For SPC/E, the self-diffusion constant DH2O is 2.25 × 10−9 m2 · s−1, while the B115 and E1.10 DH2O are respectively 0.015 × 10−9 m2 · s−1 and 0.0014 × 10−9 m2·s−1. The small diffusion constant values suggest that liquids are in a glassy-like state because of their greater intermolecular interaction strengths resulting in larger order parameters (Fig. 5). Reducing the bond angle and length gradually increases the orientational freedom of the water molecule, resulting in a greater diffusion constant and smaller H-bonding lifetime.

Fig. 6. Variation of the self-diffusion constants and hydrogen bond lifetime for various H–O–H angles and O–H length models.

To explore further the sensitivity to changes in bond length (elastic model) or bond angle (bent model), the data of the H-bonding lifetime in Fig. 6 are fitted. Interestingly, we find the slope of the elastic model (43.50 ps · Å−1) is much more than that of the bent model (0.52 ps · deg−1). This implies that the H-bonding ability is tremendously sensitive to changes of the O–H length. The intermolecular interaction strength is more sensitive to the water molecular bond length alteration than to the changes in the bond angle. As will be discussed below, the local weak interaction analysis method can further provide a visual means of observing the variation in the intermolecular interaction strength.

3.4. Mean squared displacement

From the mean squared displacement (MSD) as a function of time, Fig. 7, it can be seen that the diffusion behavior has always been preserved irrespective of the reduction in H–O–H angle or O–H length. However, the B115 displays caging behavior before attaining the diffusive regime and the E1.10 obtains a pronounced caging feature throughout the entire simulation. Chatterjee and co-workers[22] have studied the diffusive regime of changing the bond angle at a range of temperatures. They found that the diffusive behavior transfer from a liquid-like diffusive (B115) to a pronounced caging feature (B120) when the temperature is at 330 K. The difference in changes of bond angle is only 5°. Nevertheless, in Fig. 7, the much stronger sensitivity to changes in the bond length is remarkable, comparing with the bond angle. The exhibition of the caging character is caused by changing the SPC/E’s bond length 0.1 Å at room temperature.

Fig. 7. Mean squared displacement (MSD) as a function of time for various H–O–H angles and O–H length models.
3.5. Local weak interaction

Weak interaction analysis between molecules or atoms, called non-covalent interaction analysis, has been recently developed by Contreras-Garcia et al.[44–47] It is based on an analysis of the electron density distribution in molecular systems in the regions of low electron density and low gradient. This approach is often referred to as reduced density gradient (RDG) analysis[37] and provides a good way to distinguish H-bond, vdW interaction, and steric effect.

Figure 8 displays the average reduced density gradient (aRDG) and thermal fluctuation index (TFI) for various bond length and angle models. In Figs. 8(a) and 8(c), two blue ellipses near the two hydrogens act as the donors to form the ability of H-bonds with neighboring water molecules. Also a significant area of isosurface above the nearest oxygen atom serves as an indicator when the water acceptors form H-bonds with neighboring water molecules. This is shown by the colors of both regions fading with decreasing the O–H length and H–O–H angle indicating the reduction of electrostatic attraction between water molecules, and also further confirms the above mentioned suggestion, i.e., the strengths of forming an H-bond between water molecules are changed by altering the water structure. A similar argument applies to the colour changes of TFI in Figs. 8(b) and 8(d). Comparing the color variations shows that water molecules as a donor forming H-bonds with other water molecules are more stable than those as an acceptor forming H-bonds. Bartha[48] found that when H-bonds are formed between water molecules, the electronic distributions of both the donor and acceptor change. The charge distribution in the lone pair electron and hydrogen regions is reduced when water acts as an acceptor. This helps hydrogen atoms become new donors and inhibits their acceptability. It also explains that the change of the local density distribution change is observed in SDFs (Fig. 3). The regions between the first and second hydration shells of water molecules are shown in a greater isosurface behind the H-bond acceptor region. In Figs. 8(b) and 8(d) it is shown that the instability of the regions gradual decreases the steric effect with decreasing the bond length and angle. This can explain that in the second hydration shell the water fades with decreasing the bond length and angle, as seen in Fig. 2(a) and Fig. 3(a).

Fig. 8. (color online) Average reduced density gradient (RDG) and thermal fluctuation index (TFI) figures of different O–H length and H–O–H angle water models, where panels (a) and (c) correspond to aRDG, (b) and (d) correspond to TFI, respectively. Both aRDG and TFI isovalue are 0.25, using Multiwfn software default parameters. For aRDG, the color depth, in blue, represents the electrostatic interaction or H-bonding strength in the corresponding region. Similarly, the deeper red color represents the more intensive steric effect. Green indicates a low electron density region, corresponding to vdW interaction. For TFI, more blue (red) means TFI is smaller (larger), and thus the weak interaction in the corresponding region is more stable (unstable).

Furthermore, the green isosurface around a water molecule exhibits in which direction this water tends to interact with other waters by vdW interaction, as shown in Figs. 8(a) and 8(c). The green isosurface increases with the O–H length and H–O–H angle decreasing, due to the reduction in the H-bonds strength. This indicates an increase in the scope of vdW interaction between water molecules. In contrast, with the O–H length and H–O–H angle increasing, due to the increase in the H-bonds strength, the network structure is strengthened resulting in the vdW interaction isosuface becoming narrower. Additionally, the totally red vdW interaction region, as seen in Figs. 8(b) and 8(d), points to the water molecule being evidently unstable compared to H-bonding.

3.6. Ion–ion and ion–water correlations

In order to investigate the solution properties of water models with the geometry modified, ion solutions in different solvent environments are examined by adding NaCl ions. Use of the ion solution generally introduces two effects due to the ions on interactions: (i) it breaks the original solvent structures, making the solvent a little more mobile, and (ii) new hydration shells around ions reduce water’s diffusion constant.[49–52]

The RDFs of ion–ion and ion–water in the NaCl solutions of different solvents are given in Fig. 9, where, see panels (a) and (d), the probability of contact Na+–Cl− ion pairs formation is increased as the H-bond strength deviates from its natural state. The variation suggests that the cations and anions are most homogeneously distributed in SPC/E. In addition, the ion–water pair correlation functions for NaCl solution are presented in Figs. 9(b), 9(c), 9(e), and 9(f). From the peaks of their first shells, it can be seen that SPC/E water molecules are more easily coupled with NaCl than in any other solvent environment.

Fig. 9. RDFs between ion–ion and ion–water for NaCl solutions. Panels (a), (b), and (c) represent the H–O–H angle variation while panels (d), (e), and (f) represent the O–H length variation.
4. Conclusion

In this work, we have studied the liquid structure, dynamic, water molecule local weak interaction, and NaCl’s solution behavior of modified water models at 1 bar and 298.15 K. These models, based on alteration of the SPC/E model, have O–H lengths ranging from 0.90 Å to 1.10 Å (the H–O–H angles do not change) and H–O–H angles ranging from 90° to 115° (the O–H length changes from 0.82 Å to 1.08 Å). Both types of models have the same atomic charge and LJ parameters as the original SPC/E model.

To summarize, we find that the properties of the liquids undergo dramatic change caused by altering the water molecule geometry structure. Increasing the H–O–H angle and O–H length leads to: the structure of liquid becoming more highly-structured, raising the order parameter and QT, longer H-bond lifetime, and shorter atoms distance adjoining water molecules. In contrast, decreasing the H–O–H angle and O–H length results in a gradual changing of the liquid structure into being more like a Lennard-Jones liquid with spherical symmetrical local arrangement. The liquid exhibits a glassy-like state feature, when the SPC/E’s O–H length is increased by 0.1 Å. All these variations can be primarily attributed to the alteration of the intermolecular interaction strength between water molecules. However, the change in the intermolecular interaction strength is very sensitive to changes in the bond length.

Using the weak interaction and thermal fluctuation analysis methods, some interesting regions are obtained around a central water molecule. The scope of vdW (van der Waals) interaction with neighboring water molecules is reduced as the bond length and angle are decreased, and vice versa. A water molecule as a donor forms H-bonds with other water molecules leading to more stability, whereas the stability slightly reduces when the water acts as an H-bond acceptor. In addition, for the NaCl dissolving process, SPC/E water is an ideal solvent which dissolves ions and more easily couples with Na+ and Cl− ions than in other solvent environment models.

This study helps gain further fundamental understanding of the property changes of a water molecule by altering the water molecule geometry structure. Ongoing research will investigate the conformation transitions of large molecules such as DNA, protein, etc., in these modified models.

Reference
[1] Davis J G Gierszal K P Wang P Ben-Amotz D 2012 Nature 491 582
[2] Stiopkin I V Weeraman C Pieniazek P A Shalhout F Y Skinner J L Benderskii A V 2011 Nature 474 192
[3] Ji M Odelius M Gaffney K J 2010 Science 328 1003
[4] Chen B Ivanov I Klein M L Parrinello M 2003 Phys. Rev. Lett. 91 215503
[5] Vega C Sanz E Abascal J L F 2005 J. Chem. Phys. 122 114507
[6] Harrington S Poole P H Sciortino F Stanley H E 1997 J. Chem. Phys. 107 7443
[7] Horne R A Johnson D S 1966 J. Phys. Chem. 70 2182
[8] Tan M L Fischer J T Chandra A Brooks B R Ichiye T 2003 Chem. Phys. Lett. 376 646
[9] Gu B Zhang F S Wang Z P Zhou H Y 2008 Phys. Rev. Lett. 100 088104
[10] Shen X Gu B Che S A Zhang F S 2011 J. Chem. Phys. 135 034509
[11] Shen H Cheng W Zhang F S 2015 RSC Adv. 5 9627
[12] Huang J Lopes P E M Roux B MacKerell A D 2014 J. Phys. Chem. Lett. 5 3144
[13] Vedamuthu M Singh S Robinson G W 1994 J. Phys. Chem. 98 2222
[14] Yan Z Buldyrev S V Stanley H E 2008 Phys. Rev. E 78 051201
[15] Nilsson A Pettersson L G M 2015 Nat. Commun. 6 8998
[16] Lynden-Bell R M Giovambattista N Debenedetti P G Head-Gordon T Rossky P 2011 Phys. Chem. Chem. Phys. 13 2748
[17] Bergman D L Lynden-Bell R M 2001 Mol. Phys. 99 1011
[18] Lynden-Bell R M Head-Gordon T 2006 Mol. Phys. 104 3593
[19] Lynden-Bell R M Youngs T G A 2006 Mol. Simul. 32 1025
[20] Lynden-Bell R M Debenedetti P G 2005 J. Phys. Chem. B 109 6527
[21] Zhang F S Lynden-Bell R M 2005 Phys. Rev. E 71 021502
[22] Chatterjee S Debenedetti P G Stillinger F H Lynden-Bell R M 2008 J. Chem. Phys. 128 124511
[23] Gu B Zhang F S Wang Z P Zhou H Y 2008 J. Chem. Phys. 129 184505
[24] Berendsen H J C Grigera J R Straatsma T P 1987 J. Phys. Chem. 91 6269
[25] Weerasinghe S Smith P E 2003 J. Chem. Phys. 119 11342
[26] Bussi G Donadio D Parrinello M 2007 J. Chem. Phys. 126 014101
[27] Berendsen H J C Postma J P M van Gunsteren W F DiNola A D Haak J R 1984 J. Chem. Phys. 81 3684
[28] Ryckaert J P Ciccotti G Berendsen H J C 1977 J. Comput. Phys. 23 327
[29] Darden T York D Pedersen L 1993 J. Chem. Phys. 98 10089
[30] Petersen H G 1995 J. Chem. Phys. 103 3668
[31] Hess B Kutzner C van der Spoel D Lindahl E 2008 J. Chem. Theory Comput. 4 435
[32] Lu T Chen F 2012 J. Comput. Chem. 33 580
[33] Humphrey W Dalke A Schulten K 1996 J. Mol. Graphics 14 33
[34] Kumar R Schmidt J R Skinner J L 2007 J. Chem. Phys. 126 204107
[35] Kalinichev A G Bass J D 1994 Chem. Phys. Lett. 231 301
[36] Svishchev I M Kusalik P G 1993 J. Chem. Phys. 99 3049
[37] Errington J R Debenedetti P G 2001 J. Chem. Phys. 409 318
[38] Torquato S Truskett T M Debenedetti P G 2000 Phys. Rev. Lett. 84 2064
[39] Sharma R Chakraborty S N Chakravarty C 2006 J. Chem. Phys. 125 204501
[40] Truskett T M Torquato S Debenedetti P G 2000 Phys. Rev. E 62 993
[41] Abascal J L F Vega C 2007 J. Phys. Chem. 111 15811
[42] Abascal J L F Vega C 2007 Phys. Chem. Chem. Phys. 9 2775
[43] Jabes B S Nayar D Dhabal D Molinero V Chakravarty C 2012 J. Phys: Condens. Matter 24 284116
[44] Johnson E R Keinan S Snchez P M Garca J C Cohen A J Yang W 2010 J. Am. Chem. Soc. 132 6498
[45] Contreras-Garcia J Johnson E R Keinan S Chaudret R Piquemal J. P Beratan D N Yang W 2011 J. Chem. Theory Comput. 7 625
[46] Garcia J C Calatayud M Piquemal J P Recio J M 2012 Comput. Theor. Chem. 998 193
[47] Roza A O Johnson E R Garcia J C 2012 Phys. Chem. Chem. Phys. 14 12165
[48] Bartha F Kapuy O Kozmutza C Van Alsenoy C 2003 J. Mol. Struct. (Theochem) 666 117
[49] Lyubartsev A P Laaksonen A 1996 J. Phys. Chem. 100 16410
[50] Gu B Zhang F S Huang Y G Fang X 2010 Chin. Phys. B 19 036101
[51] Kropman M F Bakker H J 2001 Science 291 2118
[52] Hribar B Southall N T Vlachy V Dill K A 2002 J. Am. Chem. Soc. 41 12302